home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / sph_4pnt.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  119 lines
  1. ; $Id: sph_4pnt.pro,v 1.10 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       SPH_4PNT
  8. ;
  9. ; PURPOSE:
  10. ;       Given four 3-dimensional points, this procedure returns the
  11. ;       center and radius necessary to define the unique sphere passing
  12. ;       through those points.
  13. ;
  14. ; CATEGORY:
  15. ;       Analytic Geometry.
  16. ;
  17. ; CALLING SEQUENCE:
  18. ;       SPH_4PNT, X, Y, Z, Xc, Yc, Zc, R
  19. ;
  20. ; INPUTS:
  21. ;       X: A 4-element vector containing the X coordinates of the points.
  22. ;       Y: A 4-element vector containing the Y coordinates of the points.
  23. ;       Z: A 4-element vector containing the Z coordinates of the points.
  24. ;
  25. ;    Note: X, Y, and Z should be floating-point or double-precision
  26. ;          vectors.
  27. ;
  28. ; OUTPUTS:
  29. ;       Xc: The sphere's center x-coordinate. 
  30. ;       Yc: The sphere's center y-coordinate.
  31. ;       Zc: The sphere's center z-coordinate.
  32. ;       R:  The sphere's radius.
  33. ;
  34. ; RESTRICTIONS:
  35. ;       Points may not coincide.
  36. ;
  37. ; EXAMPLE:
  38. ;       Find the center and radius of the unique sphere passing through
  39. ;       the points: (1, 1, 0), (2, 1, 2), (1, 0, 3), (1, 0, 1)
  40. ;       
  41. ;       Define the floating-point vectors containing the x, y and z 
  42. ;       coordinates of the points. 
  43. ;         X = [1, 2, 1, 1] + 0.0
  44. ;      Y = [1, 1, 0, 0] + 0.0
  45. ;      Z = [0, 2, 3, 1] + 0.0
  46. ;
  47. ;       Compute the sphere's center and radius.
  48. ;         SPH_4PNT, X, Y, Z, Xc, Yc, Zc, R
  49. ;
  50. ;       Print the results.
  51. ;         PRINT, Xc, Yc, Zc, R
  52. ;
  53. ; MODIFICATION HISTORY:
  54. ;       Written by:  GGS, RSI, Jan 1993
  55. ;       Modified:    GGS, RSI, March 1994
  56. ;                    Rewrote documentation header.
  57. ;                    Uses the new Numerical Recipes NR_LUDCMP/NR_LUBKSB.
  58. ;       Modified:    GGS, RSI, November 1994
  59. ;                    Changed internal array from column major to row major.
  60. ;                    Changed NR_LUDCMP/NR_LUBKSB to LUDC/LUSOL
  61. ;       Modified:    GGS, RSI, June 1995
  62. ;                    Added DOUBLE keyword.
  63. ;       Modified:    GGS, RSI, April 1996
  64. ;                    Modified keyword checking and use of double precision.
  65. ;-
  66.  
  67. PRO SPH_4PNT, X, Y, Z, Xc, Yc, Zc, R, Double = Double
  68.  
  69.   ON_ERROR, 2
  70.  
  71.   if N_PARAMS() ne 7 then $
  72.     MESSAGE, "Incorrect number of arguments."
  73.  
  74.   TypeX = SIZE(X) & TypeY = SIZE(Y) & TypeZ = SIZE(Z)
  75.  
  76.   if TypeX[TypeX[0]+2] ne 4 or $
  77.      TypeY[TypeY[0]+2] ne 4 or $
  78.      TypeZ[TypeZ[0]+2] ne 4 then $
  79.      MESSAGE, "X, Y, and Z coordinates are of incompatible size."
  80.  
  81.   ;If the DOUBLE keyword is not set then the internal precision and
  82.   ;result are determined by the type of input.
  83.   if N_ELEMENTS(Double) eq 0 then $
  84.     Double = (TypeX[TypeX[0]+1] eq 5 or $
  85.               TypeY[TypeY[0]+1] eq 5 or $
  86.               TypeZ[TypeZ[0]+1] eq 5)
  87.  
  88.   if Double eq 0 then A = FLTARR(3,3) else A = DBLARR(3,3)
  89.  
  90.   ;Define the relationships between X, Y and Z as the linear system.
  91.     for k = 0, 2 do begin
  92.       A[0, k] = X[k] - X[k+1]
  93.       A[1, k] = Y[k] - Y[k+1]
  94.       A[2, k] = Z[k] - Z[k+1]
  95.     endfor
  96.  
  97.   ;Define right-side of linear system.
  98.     Q = X^2 + Y^2 + Z^2
  99.     C = 0.5 * (Q[0:2] - Q[1:3])
  100.  
  101.   ;Solve the linear system Ay = c where y = (Xc, Yc, Zc)
  102.     LUDC, A, Index, Double = Double
  103.   ;Solution y is stored in C
  104.     C = LUSOL(A, Index, C, Double = Double)
  105.  
  106.   ;The sphere's center x-coordinate.
  107.     Xc = C[0]
  108.  
  109.   ;The sphere's center y-coordinate.
  110.     Yc = C[1]
  111.  
  112.   ;The sphere's center z-coordinate.
  113.     Zc = C[2]
  114.  
  115.   ;The sphere's radius.
  116.     R = SQRT(Q[0] - 2*(X[0]*Xc + Y[0]*Yc + Z[0]*Zc) + Xc^2 + Yc^2 + Zc^2)
  117.  
  118. END
  119.